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Abstract. In this work we present a detailed computational study of structural 
and elastic properties of cubic Al x Ga y lni_ x _yN alloys in the framework of 
Keating valence force field model, for which we perform accurate parametrization 
based on state of the art DFT calculations. When analyzing structural properties, 
we focus on concentration dependence of lattice constant, as well as on the 
distribution of the nearest and the next nearest neighbour distances. Where 
possible, we compare our results with experiment and calculations performed 
within other computational schemes. We also present a detailed study of 
elastic constants for Al x Ga y Ini_ x _ y N alloy over the whole concentration range. 
Moreover, we include there accurate quadratic parametrization for the dependence 
of the alloy elastic constants on the composition. Finally, we examine the 
sensitivity of obtained results to computational procedures commonly employed 
in the Keating model for studies of alloys. 



1. Introduction 

Recently A1N, GaN and InN trigger a lot of interest in optoelectronics, mostly 
because of their electronic structure, which makes them very promising materials for 
application in blue/green and UV active devices. Emitters and detectors operating 
in these spectral range can be used in many important areas such as optical data 
storage, biosensors, multimedia etc. Nitride alloys allow for continuous tuning of 
physical quantities such as band gap, lattice parameters, mobility etc., to reach the 
suitable values for desired applications. 

Despite both theoretical and experimental efforts, many basic properties of nitride 
alloys are not yet sufficiently well understood. In this work, we focus on calculations 
of structural and elastic properties of quaternary Al x Ga v Ini_ x _ y N alloys, since they 
offer the largest possibility of tuning. In particular, we examine the morphology 
of alloys concentrating on atomic distance distribution between the nearest and the 
next nearest neighbours. They have important influence on the electronic structure 
of alloys (see e.g. [H, This knowledge is crucial for application purposes. In 
the case of simpler ternary alloys AlGaN @, AlInN @, GalnN 0, i, SS, H, 
the bond lengths distribution has been obtained by extended X-ray absorption fine 
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structure spectroscopy (EXAFS), however, for quaternary Al x Ga y In 1 _ x _ y N there are 
no available experimental data. As far as clastic properties of alloys are concerned, 
their correct description is also very important issue in modelling lasing from quantum 
wells, e.g. within k ■ p or similar continuous models. The nonlinear effects in 
elasticity of nitrides and their influence on properties of devices have recently attracted 
considerable attention [l(J EH EH E3] • However, there is not much known about the 
detailed dependence of elastic properties on composition. In this paper, we present a 
computational study of the second-order elastic constants, cy, for Al x Ga y Ini_ x _ y N 
quaternary alloys within the whole concentration range. This could provide very useful 
insight from viewpoint of device modelling. 

In the present study of quaternary nitride alloys, our main computational tool is 



valence force field approach (VFF) developed by Keating jl4j . Even though it was 
developed over forty years ago, it is still important ingredient of multiscale models, 
particularly where large number of alloy configurations needs to be handled. For 
nitrides, the Keating VFF model has been recently used to examine a plethora of 



physical phenomena, such as phonon spectra in bulk and in supcrlattices [151 . 116 



structural properties of ternary bulk alloys [18|, LHJ, |2CJ, 1211, |22| an( i their nanowires 



17]. 

231, 



stability of different all oy p has es | 22lj and also in numerous studies of thermodynamics 
of ternary 24, 2f| 2^, 27, 2|| 29{ and quaternary [3(| nitride alloys. Generally, the 
Keating VFF model is also a method of choice, where the atomic positions are needed 
as external input for electronic structure modelling. This is the case for methods 
atomistic in nature, but not entirely based on first principles, such as semiempirical 
tight-binding or empirical pseudopotential methods, that are commonly used for 
studies of low dimensional semiconductor structures. Therefore, to contribute to 
further development and validation of the model, we also pay attention to practical 
aspects of VFF usage. We compare the distribution of the nearest neighbour and next 
nearest neighbour distances resulting from Keating VFF model with those obtained 
from accurate quantum mechanical formalism, which is a good test of VFF model 
reliability. We also examine the influence of so called mixing rule used to obtain VFF 
parameters for alloys, which shows how strongly this could influence prediction of this 
model. 

The paper is organized as follows. In section [2] we briefly recall basic facts about 
Keating model and present the employed set of parameters. Section[3]provides detailed 
overview of structural properties for quaternary nitride alloys, it includes lattice 
constants as well as the distribution of the nearest neighbour and the next nearest 
neighbour distances. In this part we also compare the results of Keating VFF with 
DFT findings. In section [4] the VFF results for elastic constants of Al x Ga y Ini_ x _ y N 
alloys in the whole concentration range are presented. Section [5] deals with the 
computational procedures, so-called mixing rule used to obtain VFF parameters of 
alloys and the effect of finite supcrccll size. Finally, the paper is summarized in 
section [5] 

2. Keating valence force field model and its parametrization 

Keating [l4| , on the basis of general symmetry considerations, derived potential energy 
model for zinc-blende type crystals in the following form 

vw 2 ,...)=£ £ iH-H-4) 2 + (i) 

i jeNN(i) V 
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^2 ^2 Wd^-d ^ Vlj ' V%k ~ d%]d%k cos °o) 2 ■ 

i 3,fceJVN(i) % 3 lk 

Here, = r$ — rj , where r j , r j denote the position of i-th and j-th atoms respectively, 
dij denotes equilibrium distance between atom i and j, NN(i) represents the set of 
four nearest neighbours of i-th atom, ctij and /3ijk stand for force constants. 

In the studies applying VFF to nitrides, very often a set of parameters proposed 
by Kim et al [3l| is used. However, during the last years theoretical suggestions 
appeared, that propose refined procedure for determination of the force constants in 
VFF model [2l| . On the other hand, the overall improvement of accuracy of the first- 
principles computational methods has been also achieved. Therefore, we have decided 
to recalculate the force constants and bring them to the state-of-the-art. 

The parametrization of the Keating VFF model starts from determination of the 
force constants a and f3 for bulk zinc-blende compounds. This can be done knowing 
the elastic constants of the bulks. As shown by Keating (l4| . the a and [3 are given 
by analytic formula as a function of the elastic constants en and c\2- In the Keating 
model, the third elastic constant, i.e., C44, is related to en and C12 as follows 

2c4 4 (cn + C12) _ ^ ^ 
(en -ci 2 )(cu +3ci 2 ) 

Therefore, in the standard procedure to determine the a, and /3, C44 is not taken 
into account. However, the equation ^ for elastic constants resulting from Keating 
model is not very well satisfied for nitrides. To improve this approach, Grosse and 
Neugebauer [2lJ proposed an alternative method. They suggested to determine a 
and f3 by the least-squares fit to all three elastic constants. It turns out that such 
fitting approach ensures more uniform spreading of error and generally leads to 
better results. Therefore, we follow this approach in our work. Since the values 
of elastic constants for zinc-blende nitrides have not been measured so far, we relay 
on theoretical predictions from the DFT based calculations. In our study, we take 
arithmetic average of clastic constants obtained in the calculations within the DFT 
LDA and DFT GGA approximations to determine the force constants a and ft. The 
clastic constants obtained through the averaging over these two theoretical schemes 
should be closer in value to experimental ones. It follows from the DFT LDA and 
DFT GGA tendencies to, respectively, overestimate and underestimate the stiffness 
of a material. The values of Cij calculated within GGA formalism have been taken 
from our previous work [12j , and the Cp yalues have been calculated within DFT LDA 
approximation using VASP package 32, 33, HH and the projector augmented wave 
method |35i. The energy cutoff was set to 800 cV and the llxllxll Monkhorst-Pack 
mesh |36| has been employed for the Brillouin zone integrals. The values of elastic 
constants used for the fitting procedure arc summarized in the table [1] whereas the 
final set of employed parameters is shown in the table [2] 

Now we are in the position to determine parameters ctij and ftijk for calculation 
of Al x Ga y Ini_ x _ y N quaternary alloys. The values of the force constants ay are 
directly taken as force constants a of binary compound consisting of atomic species 
i and j. In the case of ftijk three types of atoms can be involved, and obviously the 
parametrization, carried out for binary materials cannot determine these constants. In 
such cases we employ arithmetic mixing rule, i.e., take arithmetic average of suitable 
binary parameters, as it has been applied in many previous works, e.g. flil. I25I l26j|. 
Specifically, this means that /^Ga.N.in = (/?Ga,N.Ga + An.N,in)/2, and so on. However, 
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Table 1. Elastic constants used to parametrize Keating model. 





Cll 


C12 


c 44 




A1N 


301 


167 


188 


LDA 








282 


149 


179 


GGA, 


see 


[12] 


GaN 


288 


159 


161 


LDA 






252 


129 


147 


GGA, 


sec 


[12] 


InN 


184 


127 


84 


LDA 






159 


102 


78 


GGA, 


see 


[12] 



there is also possibility to use geometric average instead (see e.g. Schabel and Martins 



37| or Saito and coworkers |20(). We investigate the role of employed mixing rule 



later on (see section [5] for details). 



Table 2. Set of parameters for Keating VFF model employed in the present 
study. 





d[A] 


a [N/m] 


P [N/m] 


A1N 


1.894 


79.91 


19.73 


GaN 


1.950 


76.25 


17.80 


InN 


2.165 


62.07 


9.68 



3. Structural properties from VFF model and their comparison with 
DFT calculations 

In this section, we analyze composition dependence of structural properties and 
geometry of Al x Ga y Ini_ x _ y N alloys (i.e., lattice constants, the nearest and next 
nearest neighbour distances). The distribution of bond lengths in an alloy can be 
also extracted from EXAFS experiments, which is a useful crosscheck for the theory. 
The local geometry of alloy is an important issue, since it influences the electronic 
structure, see e.g. recent calculations of Gorczyca et al [H 0] for nitride alloys. 
The Keating VFF model is particularly suitable for the task of establishing local 
environment of random alloys. It enables calculations with large supercells, which in 
turn guarantee reasonably good sampling and randomness of alloys. There is fairly 
long history of employing VFF models to determine alloy geometry for various types of 
materials. For example, Cai and Thorpe studied relaxation patterns in general ternary 
A\- X B X C [HI and quaternary Ai^ x B x C\- y D y alloys [39[ using Kirkwood model very 



similar to Keating VFF approach. Schabel and Martins [37| gave detailed overview 
of structure for very broad range of semiconducting alloys within Keating VFF. Their 
work, however, does not include nitrides. The structure of ternary nitride alloys has 
been also studied by many other authors 0, [H, H3| . In this work, we present the first 
VFF study of structure for quaternary nitride alloys. We compare our results for bond 
length distribution with recent calculations of Marques et al [40| carried out within 
generalized quasichemical approximation (GQCA). Where possible we also compare 
our findings with experimental data, obtained using EXAFS approach. To shed some 
light on accuracy of Keating VFF, we also compare the structural properties of this 
force field model with accurate DFT calculations. The latter treatment allows only 
for calculations with moderate sizes of supercells. On the other hand DFT description 
of interactions between atoms is otherwise very complete. Therefore, such comparison 
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can reveal potential weak points of fully local classical description. Since VFF model is 
widely used to provide geometry to semi-empirical calculations of electronic structure 
handling large systems (e.g. tight-binding or empirical pseudopotcntial schemes), the 
question of its accuracy is vivid. Obviously, the quality of such multiscalc approach 
depends on accuracy of input structures, however, to the best of our knowledge, no 
detailed comparison of VFF with more accurate models (e.g. DFT) has been so far 
presented in the literature. 



3.1. Computational details 

For VFF simulations, we used 18 x 18 x 18 zinc-blende cubic supcrcells containing 46656 
atoms with random distribution of cations. Wc optimized both the atomic positions 
and the lattice constant of the cubic cell. For the DFT simulations, we employed 



VASP package [32|, |33j, 13J] . We used local density approximation for exchange and 
correlation functional according to Ccperley and Alder [4lJ . The projector augmented 
wave method was used in its variant available in the VASP code [35[ . The calculations 
were performed for 3 x 3 x 3 zinc-blende cubic cells containing 216 atoms. Also in this 
case cation distribution among sites was random. The energy cutoff was set to 550 
eV. For Brillouin zone sampling, only Y point was used due to the large supercell size. 



3.2. Lattice constant 

In the early 20th century Vegard noticed that lattice constant for alloys can be 
calculated using linear interpolation between lattice constants of constituents (42l |. 
Nitride alloys follow this so-called Vegard's law, which has been confirmed by a series 
of experimental findings, see e.g. recent experiment for GalnN nanowires, where the 
linear dependence of lattice constant was observed in the whole concentration range 
43|. 

The results of our simulations for lattice constant are presented in figure [T] One 
can easily observe that DFT LDA formalism reproduces very well the linear Vegard- 
likc behaviour, even though the LDA approximation is well known to systematically 
underestimate the values of lattice constants. On the other hand the Keating VFF 
predicts some bowing in the dependence of lattice constant on concentration. This is in 
accordance with the work of Thorpe and coworkers [44j , who analyzed the Vegard's law 
for a simple VFF-like model. In their case they showed that Vegard's law is obtained 
when force constants disorder is neglected. Therefore, since in our calculations force 
constant disorder is included (i.e., a and (3 depend on the types of atoms considered), 
it is possible to expect deviations from linearity. This is a feature of the Keating VFF 
method, however, as can be seen from figure [U its magnitude is not very large. 



3.3. Nearest neighbours distance 

Results of calculations of the nearest neighbour distance dependence on concentration 
of alloy constituents for Al x Ga y Ini_ x _ y N are presented in figure [2] It is well known 
that even though lattice constants obey Vegard's law, the individual bond lengths 
do not follow this simple rule. The dependence of bond lengths on concentration is 
usually also linear, however, bond lengths remain much closer to their original bond 
length in bulk binary material, rather then to the average bond length predicted by 
Vegard-like law (see e.g. 4a, |46j). This linear dependence of the bond lengths in the 
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Figure 1. Comparison of lattice constants dependence on the alloy composition 
for Al x Ga y Ini_ x _ y N obtained using Keating VFF (left panel) and DFT LDA 
(right panel). Points correspond to results of calculations, solid lines are only 
to guide the eye. Dashed lines denote prediction of Vegard's law. Note that for 
convenient comparison the scale on both graphs is the same. 



alloy on the concentration of Al and Ga cations, x and y respectively, can be described 
by the linear form 

d(x, y)=d + Ax + By. (3) 

The coefficients A and B in equation ([3j have been determined by fitting the 
bond lengths obtained from both Keating and LDA schemes. These linear fits are 
represented by dashed lines in figure [21 The small values of coefficients A and B 
(see table [3] gathering results of fitting procedure for all types of bonds) allow even 
for rough approximation of d(x,y) by do, just confirming the statement about bond 
lengths made above. 

On the basis of the presented plots (figure[2|) for bond lengths in Al x Ga y Ini_ x _ y N 
alloy, one can compare findings of Keating VFF model with DFT LDA approach. 
Results of both computational methods can be reasonably well described by linear 
dependence of bond lengths on cation concentrations given by equation ([3]). Naturally, 
since the DFT LDA calculations have been performed for much smaller cell sizes 
(histograms were generated on the basis of two 216 atoms supercell calculations) the 
data have larger statistical error than the Keating VFF computations. This effect 
is particularly pronounced in these parts of the graphs where low concentration of 
considered cation is present in the sample. Second evident observation is that DFT 
LDA bond lengths are systematically shifted towards the lower values than Keating 
VFF. This is related to well known LDA flaw to underestimate the lattice constants. 
It is also worth noticing that results of previous computations by Marques et al (4p| . 
which have been carried out within generalized quasichemical approximation (GQCA), 
are in particularly good agreement with findings of Keating VFF model (see table [3] 
for details). 

Another interesting aspect is to compare the probability distribution profiles 
predicted by both DFT LDA and Keating VFF models. Two sample histograms 
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generated by both models are depicted in figure[3J One can see that both schemes lead 
to results that are in satisfactory agreement with each other. It is again visible that 
DFT LDA systematically underestimates the bond lengths, so the peaks corresponding 
to each bond type are shifted towards lower values. One can also notice that 
Keating VFF model predicts larger peak widths than the DFT LDA approach does. 
Particularly interesting are cases of high In concentration. There, one can observe 
rather large lattice deformations caused mostly by the considerable lattice mismatch 
between InN and GaN or A1N. This can be considered as a test of the Keating VFF 
which has been parametrized on the basis of elastic constants aj, i.e., taking into 
account effects of second-order in deformation strain, which, in principle, describe the 
material behaviour in the regime of small deformations. However, also for these cases, 
the agreement is reasonable, in spite of neglecting higher order contributions to the 
elastic energy (see e.g. results for Alo.17Gao.17Ino.66N presented in figured]). 

Table 3. Results of the linear fits to average nearest neighbours distances for 
Al x Ga y Ini_ x _ y N quaternary alloys. 



Bond 


Model 


Fit Function 


Error [%] 


Ref. 


Al-N 


VFF 
DFT LDA 
GQCA 


1.9496 - 0.05496 x - 0.03893 y 
1.9279 - 0.04419 x - 0.04314 y 
1.9435 - 0.05090 x - 0.05455 y 


0.10 
0.26 
0.19 


This work 
This work 
[40] 


Ga-N 


VFF 
DFT LDA 

GQCA 


1.9959 - 0.06322 x - 0.04573 y 
1.9762 - 0.04367 x - 0.04470 y 
1.9919 - 0.05985 x - 0.06050 y 


0.10 
0.27 
0.35 


This work 
This work 
[40] 


In-N 


VFF 
DFT LDA 

GQCA 


2.1676 - 0.09899 x - 0.07485 y 
2.1461 - 0.05276 x - 0.05088 y 
2.1573 - 0.08292 x - 0.07570 y 


0.12 
0.21 
0.55 


This work 
This work 
[40] 



Finally, it is interesting to examine, how the presented theoretical results agree 
with experiments. We are not aware of any experimental results for bond length 
distribution in Al x Ga y Ini_ x _ y N alloys, however, there is a significant body of 
experimental data for ternaries Al x Gai_ x N, Al x Ini_ x N, and Ga x Ini_ x N. Comparison 
of the experimental bond lengths in ternary alloys with predictions of the Keating VFF 
model is presented in figure |H Even though measurements are performed for wurtzitc 
samples, the comparison with the theoretical predictions for cubic phases is reasonable, 
owing to the local similarity of both crystallographic phases. In wurtzite structures, 
the nearest neighbours bond elongations in the direction of crystallographic c-axis are 
typically of the order of 0.01 A [3, EtJ- This is very similar to a typical experimental 
error bars or bond length spread resulting from disorder in alloy (compare histogram 
in figure [3|) . As can be seen in figure IU the comparison of available data with our 
theoretical results reveals good agreement. It could be, of course, interesting to see 
how this relates to the structural properties of cubic Al x Ga y Iiii_ x _ y N which growth 
using molecular beam epitaxy was recently reported (48j . For completeness, in figure 
[4]we have also included the DFT findings which again show the well known systematic 
tendency to underestimate the bond lengths. 
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Figure 2. The dependence of the average bond length on alloy composition 
for Al x Ga y Ini_ x _ y N obtained using Keating VFF and DFT LDA schemes. 
Points correspond to results of calculations, solid lines are only to guide the eye. 
Dashed lines denote linear fits presented in the table [3] Note that for convenient 
comparison the scale on Keating VFF and DFT LDA graphs is the same. 
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Figure 3. The nearest neighbour distance distribution resulting from Keating 
VFF and DFT LDA for Alo.50Gao.33Ino.17N and Alo.17Gao.17Ino.66N random 
alloys. Vertical dashed lines denote nearest neighbour distance in pure A1N, GaN 
and InN respectively as predicted by both presented models. 
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Figure 4. Comparison of our theoretical results obtained within both DFT LDA 
and Keating VFF models with various experimental data. There are used data 
from the followingexperimental: for Al x Gai_ x N Expt. 1 - Miyano et al Q, 
Expt. 2 - Yu et al [J, for Al x Im_ x N Expt. 3 - Katsikini et al [J, for Ga x Im_ x N 
Expt. 4 - MBE samples of Kachkanov et al [1, Expt. 5 - MOCVD samples of 
Kachkanov et al Expt. 6 - Katsikini et a/Q, Expt. 7 - Katsikini et al ff|. 
The results of previous theoretical predictions within GQCA [4o| are presented 
for completeness. 
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3.4- Next nearest neighbour distance 

Here we analyse the closest distances between the atoms of the same (cationic or 
anionic) sublattice, i.e., the next nearest neighbour distance of an anion or cation. 
The results obtained within both VFF and DFT models for a nitrogen-nitrogen pairs 
are presented in figure[5l Qualitatively, the results for the smallest distances of cationic 
pairs are very similar, exhibiting dominantly linear dependence on concentration of 
constituents. The coefficients A and B that determine this linear dependence (see 
equation ([3])) arc presented in the tabled for both cation-cation and nitrogen-nitrogen 
average distance. The maximum error of the fit is also included there. Since in modern 
EXAFS experiments it is possible to measure the cation-cation distances for every 
distinct pair separately, we have also performed the fits for all possible combinations 
of cation-cation distances, namely Al-Al, Al-Ga, Al-In, Ga-Ga, Ga-In, In-In (see table 
[5]). Again, in all cases linear model provides satisfactory description. 



Table 4. Results of linear fit to the average next nearest neighbour cation- 
cation (CT - CT) and nitrogen-nitrogen (N - N) distances for Al x Ga y ini_ x _yN 
quaternary alloys. 



Distance 


Model 


Fit Function 


Max. Error [%] 


CT-CT 


VFF 
DFT LDA 


3.5147 
3.4934 


- 0.43234 x - 0.33500?/ 

- 0.42243 x - 0.34107 y 


0.59 
0.18 


N-N 


VFF 
DFT LDA 


3.5200 
3.4991 


- 0.43517 x - 0.33931?/ 

- 0.42487 x - 0.34562?/ 


0.44 
0.07 



Table 5. Results of linear fits to the average next nearest neighbour distances 
between various pairs of cations for Al x Ga y Ini_ x _ y N quaternary alloys. 



Distance 


Model 


Fit Function 


Max. Error [%] 


Al-Al 


VFF 
DFT LDA 


3.3660 - 0.27416 x - 0.20666 y 
3.3663 - 0.29078 x - 0.23771 y 


0.15 
0.78 


Al-Ga 


VFF 
DFT LDA 


3.3865 - 0.28385 x - 0.21459 y 
3.3798 - 0.29593 x - 0.23651 y 


0.07 
0.24 


Al-In 


VFF 
DFT LDA 


3.4599 - 0.32354 x - 0.24685 y 
3.4286 - 0.30625 x - 0.23513 y 


0.05 
0.31 


Ga-Ga 


VFF 
DFT LDA 


3.4057 - 0.29338 x - 0.22179 y 
3.3922 - 0.29127 x - 0.23931 y 


0.18 
0.43 


Ga-In 


VFF 
DFT LDA 


3.4771 - 0.33479 x - 0.25554 y 
3.4474 - 0.30469 x - 0.24991 y 


0.06 
0.26 


In-In 


VFF 
DFT LDA 


3.5362 - 0.38725 x - 0.29734 y 
3.4983 - 0.31251 x - 0.26031 y 


0.11 
0.53 



In addition to the dependence of average cation-cation and anion-anion distance 
on concentration, we have also analyzed the shape of distribution. Sample histograms 
are depicted in figure El When inspecting presented graphs, one can notice that 
cationic and nitrogen sublattices relax in a different manner. The behaviour of cation- 
cation distribution is similar to virtual crystal exhibiting unimodal shape. At the 
same time nitrogen-nitrogen distance is much more distorted with respect to single- 
peaked virtual crystal picture. The shape of the distribution is multimodal. The 
peaks correspond to three possible combinations of N-N pair joined by Al, Ga or In 
atom respectively. Since the equilibrium distances of Al-N and Ga-N are very similar 
the pairs of N-N joined by Al and Ga form common maximum on the presented 
plot. These findings agree very well with the same behaviour reported for ternaries 
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Comparing bond length histograms obtained within the VFF and DFT 
schemes (see figure , it is interesting to note that the agreement in the predicted 
next nearest neighbour distances is in this case even better than for the nearest 
neighbours histograms. To some extent it is caused by the fact that the larger 
peak widths than in the case of the nearest neighbour distance distributions cause 
that systematic differences between bond length predicted by DFT and VFF are less 
pronounced. Since the VFF model includes explicit interactions only between the 
nearest neighbours, one could think that the description of second coordination shell 
would be less accurate. However, it turns out that the VFF model provides also very 
reliable predictions for cation-cation and anion-anion distances (see figure [6]). 

Even though the experimental data for quaternaries are unavailable, similarly like 
in the nearest neighbour case, there is a considerable number of results for ternaries. 
The comparison of experimental findings with our theoretical results are presented 
in figure [7] Generally, the agreement is very good, however, a few things are worth 
pointing out here. One can notice that quite often the results measured by different 
groups exhibit significant spread and, in addition, some of the experimental data 
exhibit large error bars. This underlines the fact that such measurements are on the 
verge of available experimental technique. The largest discrepancies between theory 
and experiment we observe for Al x Ini_ x N alloy (see the middle column of the graph 
[7]), however, there is only one recent report [5[ known to us, which deals with the 
structure of this material. To elucidate the matter, more experimental data would be 
very helpful. One has also to bear in mind that our calculation assume random cation 
distribution in the sample. If some kind of clustering for particular type of cations 
occurs, it could lead to modification of the presented results. 
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Figure 5. The nitrogen-nitrogen average next nearest neighbour distances as a 
function of composition for Al x Ga y Ini_ x _ y N obtained using Keating VFF and 
DFT LDA, for the case of nitrogen-nitrogen. Points correspond to results of 
calculations, solid lines are only to guide the eye. Dashed lines denote linear fits 
presented in the table [4] Note that for convenient comparison the scale on all 
graphs is the same. 
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Cation-cation distance for Al 17 Ga 17 ln 66 N 



N-N distance for Al 17 Ga 17 ln 66 N 
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Figure 6. Comparison of the next nearest neighbour distance distribution 
resulting from the Keating VFF and DFT LDA for Alo.17Gao.17Ino.66N. Vertical 
dashed lines denote next nearest neighbour distance in pure A1N, GaN and InN 
respectively as predicted by both presented models. 
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Figure 7. Comparison of our theoretical results for average next nearest 
neighbours distances with various experimental findings: for Al x Gai_ x N Expt. 1 
- Miyano et al H, Expt. 2 - Yu et al H, for Al x Ini_ x N Expt. 3 - Katsikini et 
al for Ga x Ini_ x N Expt. 4 - MBE samples of Kachkanov et al M, Expt. 5 - 
MOCVD samples of Kachkanov et al Expt. 6 - Katsikini et al Q, Expt. 7 - 
Katsikini et al All distances are in Angstroms. 
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4. Elastic constants for alloys calculated using Keating model 

In this section wc present results of calculations of clastic constants Cij for random 
Al x Ga y Ini_ x _ y N alloys. Calculations were carried out using Keating VFF and cover 
the whole concentration range. To extract the values of the alloy elastic constants, we 
have applied three types of strains to every alloy supercell: 

e A = [e,0,0,0,0,0], 

e B = [e, -e,0,0,0,0], (4) 
e c = [0, 0,0, 0,0, e]- 

For each type of the deformation, e was varied within the range of values 
{ — 1.0%, —0.5%, 0.5%, 1%} and the elastic energy has been calculated. Then, on the 
basis of strain energy relation 

1 6 

E = T7 C *J £ * e i' ( 5 ) 

three elastic constants en, C12 and C44 have been determined from parabolic fits to 
the energy for deformations ea, £b, ec • 

The results are presented in figure [5] We also included there the prediction of 
Vcgard-like law for elastic constants: 

cj^ d (x, y) = x cf N + y c.g aN + (l-x-y) cjf . (6) 

Exact functional forms of this equation for en, C12 and C44 are explicitly given in the 
table [6] However, after a brief analysis of graph one notices that Keating VFF 
results are not very well described by the above Vegard's law. This is particularly 
pronounced for elastic constants en and C44. The deviation for C12 is very weak, but 
this elastic constant has also the lowest discrepancy between materials A1N, GaN and 
InN. To fully describe the dependence of on composition one has to include bowing 
term Acy(x, y), which is defined as follows: 

Cij = c^ cgard (a;, y) + A dj(x, y), (7) 

Acij(x,y) = Px(l -x) + Qxy + Ry(l - y). 

After including this additional function and performing fitting procedure for P, Q 
and R, the VFF results are reproduced with accuracy much better than 1 GPa. The 
coefficients P, Q, and R for cubic elastic constants are provided in the table [S] 

The literature indicates that indeed bowing in the alloy elastic constants Cij should 
be expected. Chen and Sher in their book J4(| point out that the value of Ac^ should 
be always negative, which is the case in our studies. They argue that since the elastic 
properties of semiconductors correlate with inverse power of lattice constant and lattice 
constants for alloys follows Vegard's law, then the Cij dependence on composition 
should be sublinear. They also perform simple analysis within the framework of 
Keating VFF model showing that for simple ordered structures sublinear bowing in 
bulk modulus is present. Also our preliminary calculations within virtual crystal DFT 
pseudopotential scheme (VCA DFT, sometimes referenced as computational alchemy) 
predict a presence of the bowing term in Cij{x,y) (49[. However, one has to bear in 
mind that VCA DFT works best for alloys, when lattice mismatch of constituents is 
low. Larger mismatch, (as in the case of InN with both A1N and GaN) can introduce 
considerable inaccuracy to this model. 
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To summarize this section, we have shown that Keating VFF model predicts 
quadratic dependence of the elastic constants aj on Al x Ga y Ini_ x _ y N alloy's 
composition. This effect is in agreement with the previous literature studies [itl l49j. 
We believe that, even thought the effect is not very large, the awareness of it 
could improve description and modelling of devices based on low dimensional nitride 
structures, such as quantum dots and quantum wells. When the elastic constants 
of more common wurtzite phase of Al x GayIn-| _ Y _ V N are needed, one can obtain 
desired dependencies using Martin transformation |50j to the data gathered in tabic 
O It would be also interesting to compare presented results with other modelling 
approaches, since Keating VFF is a simple tool and does not capture many effects. 
The accurate experimental studies would be also of a great value here. 
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Table 6. Concentration dependence of the elastic constants ey for 
Al x Ga y Ini_ x _yN alloys. Accuracy of linear Vegard-like model is compared with 
fits including additional bowing term Ac;j . All data in GPa. 





Result of fit for c;j 


Max. 


difference: 
[GPa] 


Cll 


Vegard's law 

182.23 + 135.79 x + 105.70 y 








^Z.U /O) 




Vegard's law + bowing term Acn 
Acn = -24.51 x(l -x) + 39.70 xy 


- 15.44 1/(1 - 


y) 


0.3 


(0.1%) 


C12 


Vegard's law 

104.78 + 32.80 x + 25.02 y 






0.7 


(0.6%) 




Vegard's law + bowing term Aci2 
Aci2 = -2.22 x(l - x) + 3.10 - 


1.00 2/(1 — y) 




0.1 


(0.1%) 


C44 


Vegard's law 

67.00 + 77.71 x + 61.20i/ 






5.3 


(5%) 




Vegard's law + bowing term AC44 
Ac 4 4 = -21.06 x(l - x) + 33.99x1/ 


- 13.68i/(l - 


y) 


0.1 


(0.1%) 



5. Influence of mixing rule and finite supercell size 

In this section wc give a brief overview of two more technical aspects of presented 
calculations. First, we examine to what extent the results of our VFF alloy simulations 
are sensitive to the selection of mixing rule. Second, we analyze how effects of finite 
supercell size influence the structural properties. 

As already mentioned in section [2J following e.g. ljl, 25, (2(| and many other 



works, we have used arithmetic mixing rule to interpolate between three-body (3 
constants of base materials, i.e., 

a PAl.N.Al + /?Ga.N,Ga ,„s 
PAl,N,Ga — y \ ' 

and so on. Obviously this is not the only option. Therefore, important question arises 
- to what extent this choice influences the results presented so far? In order to check 
this, we carried out a set of simulations for geometric mixing, used e.g. by Schabcl 



and Martins [37[ or Saito et al [20(. This means that instead of equation (JSJ wc have: 

PAl.N.Ga = V PA1,N,A1 /?Ga,N,Ga (9) 

and so on. Other Keating VFF parameters are left unchanged. 

Sample results of this numerical experiment are presented in figure [9] Generally, 
it turns out that average lengths remain virtually unaffected by the type of mixing. 
On the left panel of figure [9j where average In-N distance is presented, one can see 
that curves for arithmetic and geometric mixing almost cover each other. This is 
because the difference there does not exceed 0.05 %. The behaviour is similar for both 
average nearest neighbours and next nearest neighbours distances. Small differences 
could be observed in the results for elastic constants. Maximum deviations were: 
Acn = 1.4 GPa (0.6%), Ac 12 = 0.5 GPa (0.4%), Ac 44 = 1.4 GPa (1.3%). It turns 
out that en and C44 are always larger in arithmetic than in geometric mixing model, 
whereas for C12, it is the other way round. Graphical comparison of C44 dependence 
on concentration in both approaches is depicted in the right panel of figure |H Even 
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though results for elastic properties depend on the type of mixing used, differences 
are really very small, strongly suggesting that both ways of mixing constants /3 lead 
to equally valuable predictions. Detailed verification which rule works better would 
be rather difficult, requiring very accurate experimental investigation. 



In-N distance for Al x Ga y ln 1xy N Elastic constant c 44 for A^Ga^n^^yN 




Figure 9. Comparison of calculations results using VFF model with different 
mixing rules - arithmetic and geometric. Left panel presents di n _it(x,y), the 
results for arithmetic and geometric mixing are so close to each other that it is 
impossible to distinguish them in the figure. The right panel presents similar 
comparison for the elastic constant c^(x,y). 

Another interesting technical aspect of presented simulations is verification to 
what extent the finite cell size influences obtained results. Particularly, one may 
wonder if results obtained on the basis of two 3x3x3 supercells per concentration in 
DFT LDA calculations (section [3]) reproduce sufficiently well the properties of random 
alloy. To verify this, we computed geometries of the same small cells used for DFT 
LDA calculations by means of VFF model. The comparison of data obtained for this 
small cells with our VFF results for 18 x 18 x 18 (46656 atoms) is presented in figure 
1101 One can see that on the sample diagram of nearest neighbour distance c^ai-n, 
even though differences are present, they have the form of typical statistical noise 
preserving the same linear trends as it has been found for large supercells. The same 
behaviour was observed for other pairs of the nearest neighbours. The analysis of the 
next nearest neighbour distances reveals even better agreement. This is because the 
second coordination shell in zinc-blende (or ideal wurtzite) contains 12 atoms, whereas 
the first consists of only 4 nearest neighbours. This increases the statistics and leads 
to smaller error for the next neighbours. On the basis of presented comparison, one 
can conclude that 3x3x3 supercells reproduce correctly the trends observed in the 
structural properties for much larger systems. 

6. Summary 

In this work we have presented computational study of structural and elastic properties 
of zinc-blende quaternary Al x Ga y Ini_ x _ y N alloys over the whole concentration range. 
Our main computational tool was Keating VFF model. We have started with 
presenting new parametrization of this model based on state-of-the- art quantum- 
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Figure 10. Comparison of calculations results using VFF model with different 
cell sizes. Small cell size corresponds to data gathered from two different 3x3x3 
cubic cells (216 atoms each), large cell corresponds to 18 X 18 X 18 cubic cell 
(46656 atoms). Right panel compares nearest neighbour distance d\i_ n(x, y) 
and left panel presents next nearest neighbour distance Dpf_^f(x,y). The latter 
results for both considered sizes are so close that they are impossible to distinguish 
in the picture. 



mechanical calculations within DFT formalism. Then we have shown the VFF results 
for lattice constant and distributions of the nearest neighbours and the next nearest 
neighbour distances. We have compared these predictions with accurate DFT LDA 
calculations for supercells of moderate size. It has turned out that the agreement 
is reasonable, which shows that simple nearest-neighbour interaction approximation 
made in Keating VFF sufficiently well captures the most important aspects of more 
accurate DFT picture. Then we have also used VFF model to examine the elastic 
constants, concluding that the composition dependence of exhibits deviation from 
Vegard-like model in form of sublinear bowing. This is in accordance with suggestions 
already made in the literature 46, We have also presented accurate quadratic 
function fits, which very well approximate the dependence of on composition 
including aforementioned bowing effect. This could be used to improve continuous 
models of nanostructurcs. Finally, we have examined the influence of mixing rules on 
VFF results. It has turned out that structural properties remain virtually unaffected 
when one uses geometric mixing instead of arithmetic one. The effect on clastic 
constants is larger, however still much lower than typical experimental error. 
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